#ifndef HEUMOESER_HPP
#define HEUMOESER_HPP

#include "edgeTraversal.hpp"
#include "sweSolver.hpp"
#include "solverBase.hpp"
#include "meshTypes.hpp"
#include "meshFunctions.hpp"
#include "magic_enum.hpp"
#include "vtkFunctions_export.hpp"
#include "esriDEMReader.hpp"

#include "solverSettings.hpp"

#include <filesystem>
namespace fs = std::filesystem;

using namespace hms;

struct heumoeser {

	hms::Time time;
	hms::Traversal traversal;
	hms::EdgeTraversal edgeTraversal;
	hms::CellTraversal cellTraversal;
	Index blocksize;
	bool
		useMuscl,
		useExport,
		useAtomicWrite;
	scalar
		maxCourant,
		stricklerCoeff;
	fs::path exportPath;

	heumoeser( hms::Time&& t, hms::SolverSettings&& s, scalar strickler ) :
		time { std::move(t) },
		traversal { s.traversal },
		edgeTraversal { s.edgeTraversal },
		cellTraversal { s.cellTraversal },
		blocksize { s.blocksize },
		// useFriction { s.useFriction },
		// useSourceFields { s.useSourceFields },
		useMuscl { s.useMuscl },
		useExport { s.useExport },
		useAtomicWrite { s.useAtomicWrite },
		maxCourant { s.maxCourant },
		stricklerCoeff { strickler },
		exportPath { std::move(s.exportPath) }
	{}

	template<typename MeshType>
	void operator()(MeshType& mesh){
		/**************************
		***** case properties *****
		***************************/

		bool
			useFriction          {true},
			useSourceFields      {true};
		scalar 
			sourceEndTime        {5.0 * 60},
			manningCoeff         {1.0/stricklerCoeff};
		LimiterFunction lf       {LimiterFunction::minMod};

		/**************************
		*********** DEM ***********
		***************************/
		scalar xOrigin     {103800.0};
		scalar yOrigin     {5244675.0};
		scalar nodatavalue {-999.0};
		
		fs::path demPath {"/home/finnamann/hms-pp/hmspp/example_case/validation_cases/heumoeser/dem/DEM-Ebnit-cellsize01m.asc"};
		DEMReader dem(
			demPath,
			xOrigin,
			yOrigin,
			nodatavalue
		);

		
		Index
			nx { mesh.nCellsX() },
			ny { mesh.nCellsY() },
			ng { mesh.nGhostLayers() };

		/**************************
		*********** mesh **********
		***************************/
		mesh.blocksize(blocksize);

		/**************************
		*********** Time **********
		***************************/
		// Time time;

		// time.startAt(0.0);
		// time.endAt(simulationTime);
		// time.maxTimeStep(maxTimeStep);
		// time.writeAtTimeIntervals(useWriteTimeInterval);
		// if ( useWriteTimeInterval ) time.writeTimeInterval(writeInterval);
		// else time.writeStepInterval(static_cast<hms::Index>(writeInterval));

		/**************************
		**** solver properties ****
		***************************/
		SWESolver<MeshType> solver;

		solver.useAtomicWrite(useAtomicWrite);
		solver.time(time);
		solver.traversal(traversal);
		solver.edgeTraversal(edgeTraversal);
		solver.cellTraversal(cellTraversal);
		solver.useFriction(useFriction);
		solver.useSourceFields(useSourceFields);
		solver.useMuscl(useMuscl);
		solver.limiterFunction(lf);
		solver.maxCourant(maxCourant);

		/**************************
		******** fields ***********
		***************************/
		solver.v() = Field{mesh, 2};
		solver.d() = Field{mesh, 1};
		solver.z() = Field{mesh, 1};
		solver.v().values() = 0.0;
		solver.d().values() = 0.0;
		solver.n() = ConstantField<Uniform, MeshType>{ mesh,Array1Ns<1>{manningCoeff} };
		solver.i() = SourceField<Uniform, MeshType>{ mesh, ArrayNs<1>{0} };
		
		// create empty rain array for complete mesh.
		ArrayXXs rain (1, (nx+2*ng)*(ny+2*ng));
		rain = 0.0;

		// assign bottom elevation from the DEM and create rainfall
		mesh.forEachDomainCell( [&](const auto& cell){
			// bottom elevation
			scalar x {cell.centroid()[0]};
			scalar y {cell.centroid()[1]};

			scalar z {dem.valueAt(x, y)};

			solver.z().values()(0, cell.index()) = z;

			// rain source
			scalar mass    { 0.0 };
			scalar rainMax { 10.0 * 0.02778e-3 };

			if( z < 1050.0 ){
				scalar ratio = (z-880.0) / (1200.0-880.0);
				mass  = interpolate(0.0, 0.33*rainMax, ratio);
			} else{
				scalar ratio = (z-1200.0) / (1600.0-1200.0);
				mass  = interpolate(0.33*rainMax, rainMax, ratio);
			}

			// fill rain array
			rain(0, cell.index()) = mass; 
		});

		ArrayXXs rainDomain (1, nx*ny);
		rainDomain = rain.reshaped<Eigen::RowMajor>(ny+2*ng, nx+2*ng).block(ng,ng,ny,nx).reshaped<Eigen::RowMajor>().transpose();
		
		// assign source field.
		solver.r() = SourceField<NonUniform, MeshType> { mesh, std::move(rainDomain) };

		/**************************
		**** boundary patches *****
		***************************/
		const BoundaryPatch& left   { mesh.boundaryPatch("left") };
		const BoundaryPatch& right  { mesh.boundaryPatch("right") };
		const BoundaryPatch& bottom { mesh.boundaryPatch("bottom") };
		const BoundaryPatch& top    { mesh.boundaryPatch("top") };

		/**************************
		*** boundary conditions ***
		***************************/
		Array1Xs bottomElevationLeft   { solver.z().domainValues(left) };
		Array1Xs bottomElevationRight  { solver.z().domainValues(right) };
		Array1Xs bottomElevationTop    { solver.z().domainValues(top) };
		Array1Xs bottomElevationBottom { solver.z().domainValues(bottom) };

		solver.z().setBoundaryCondition(left,   StaticValue{}, bottomElevationLeft,   false);
		solver.z().setBoundaryCondition(right,  StaticValue{}, bottomElevationRight,  false);
		solver.z().setBoundaryCondition(top,    StaticValue{}, bottomElevationTop,    false);
		solver.z().setBoundaryCondition(bottom, StaticValue{}, bottomElevationBottom, false);

		solver.d().setBoundaryCondition( left,   FreeOutflow_waterdepth{} );
		solver.d().setBoundaryCondition( right,  FreeOutflow_waterdepth{} );
		solver.d().setBoundaryCondition( top,    FreeOutflow_waterdepth{} );
		solver.d().setBoundaryCondition( bottom, FreeOutflow_waterdepth{} );

		solver.v().setBoundaryCondition( left,   FreeOutflow_velocity{} );
		solver.v().setBoundaryCondition( right,  FreeOutflow_velocity{} );
		solver.v().setBoundaryCondition( top,    FreeOutflow_velocity{} );
		solver.v().setBoundaryCondition( bottom, FreeOutflow_velocity{} );

		/**************************
		***** initialisation ******
		***************************/
		if(traversal == hms::Traversal::edges) {
			std::cout << "Edge traversal: "
				<< magic_enum::enum_name( solver.edgeTraversal() );
		} else {
			std::cout << "Cell traversal: "
				<< magic_enum::enum_name( solver.cellTraversal() );
		}
		
		solver.initialise(); 
		mesh.printInfo();

		Timer t;
		std::cout << std::setprecision(5) << std::fixed;
		std::cout << "Time: " << time.current() << "s";

		bool rainHasStopped {false};
		while ( !time.atEnd()){

			// stop the rain, when the time is right. maybe this does not work. if so, try values()-function of r(). 
			if( time.current() > sourceEndTime && !rainHasStopped){
				solver.r() = SourceField<Uniform, MeshType>{ mesh, ArrayNs<1>{0} };
				solver.useSourceFields(false);
				rainHasStopped = true;
			}

			// vtk export
			if( useExport && time.atWriteTime() )
				vtk::toVtk(exportPath,mesh,solver,time.nWritten());
			// evolution in time.
			solver.advance( time );
		}

		// last vtk export
		if( useExport )
			vtk::toVtk(exportPath,mesh,solver,time.nWritten());

		double elapsed { t.elapsed() };
		
		// Timing
		double execTimeSum {0};
		for ( auto execTime : solver.execTimings() )
			execTimeSum += execTime;

		std::cout
			<< "\nComputation time: " << elapsed << "s"
			<< "\nSum of execution times logged by solver: " << execTimeSum << "s"
			<< "\n\t" << 100 * solver.execTimings()[0] / execTimeSum
				<< "% for time step calculation"
				<< " (" << solver.execTimings()[0] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[1] / execTimeSum
				<< "% for fluxes (including slope)"
				<< " (" << solver.execTimings()[1] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[2] / execTimeSum
				<< "% for sources " << ( useFriction ?
					std::string{ "(including friction)" } :
					std::string{ "(friction disabled)" } )
				<< " (" << solver.execTimings()[2] << "s)"	
			<< "\n\t" << 100 * solver.execTimings()[3] / execTimeSum
				<< "% for new state"
				<< " (" << solver.execTimings()[3] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[4] / execTimeSum
				<< "% for updating boundaries"
				<< " (" << solver.execTimings()[4] << "s)"
			<< "\n";
	}
};


#endif